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Abstract 


This paper examines the accuracy and calculation speed for the magnetic field computation in an 
axisymmetric electromagnet. Different numerical techniques, based on an adaptive nonuniform grid, 
high order finite difference approximations and semi-analitical calculation of boundary conditions are 
considered. These techniques are being applied to the modeling of the Variable Specific Impulse 
Magnetoplasma Rocket. For high-accuracy calculations, a fourth-order scheme offers dramatic 
advantages over a second-order scheme. For complex physical configurations of interest in plasma 
propulsion, a second-order scheme with nonuniform mesh gives the best results. Also, the relative 
advantages of various methods are described when the speed of computation is an important 
consideration. 

1. Introduction 

Computing the axisymmetric magnetic field generated by the current in a coil is a numerical problem 
frequently encountered in the mathematical modeling of plasma flows, charged particle beams [1,2, 3], 
and other areas of physics. The problem can be approached in a variety of ways involving numerical 
integration methods as well as finite difference, finite volume, finite element, or other numerical 
methods for solving partial differential equations. All methods produce results, which may trade off 
computational speed for accuracy. The results depend very much on the order of approximation, the 
use of uniform or nonuniform grids, and other factors. For steady state solutions, the numerical 
approach may be different than for those solutions involving time dependence. 

One important area of interest involves the study of plasma behavior in the Variable Specific Impulse 
Magnetoplasma Rocket (VASIMR) engine [4], This new space propulsion concept is being studied at 
the Advanced Space Propulsion Laboratory (ASPL) of the NASA Johnson Space Center in Houston, 
Texas. Numerical simulation data are a crucial element of this research. They allow us to accurately 
predict system performance before developing costly experimental test equipment. In addition, 
numerical modeling is helpful in the design of plasma heating schemes and understanding plasma 
behavior in the VASIMR. 

A number of factors direct requirements on code accuracy and computational speed. For example, the 
measurement accuracy of magnetic field with most available laboratory instrumentation is about 0 . 1 % 
[5, 6]. Furthermore, magnetic coils can only be positioned with a limited degree of precision. These 
experimental realities limit the relative error, which can be experimentally verified to about 0 . 001 . On 
the other hand, simulating plasma instabilities and other more complex dynamics may require a high- 
accuracy numerical solution. Such accurate solutions may also be useful in developing plasma control 
algorithms for an operational device. Still another important issue pertains to solver speed, as the 
complex iterative mathematical modeling of a magnetoplasma thruster requires multiple fast 
recalculation of the magnetostatic problem. 

Finite difference and finite volume methods are widely used for solving the magnetostatic problem [7- 
9], In previous publications, researchers have approached the problem in Cartesian coordinates. We, 
in turn, use a cylindrical coordinate formulation, assuming axial symmetry, and explore a high-order 
finite volume approximation. In addition, we analyze and describe the effects of uniform vs. 
nonuniform grids with respect to accuracy, computational speed and robustness. The results are 
presented in the context of the VASIMR described briefly in the next section. 
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2. VASIMR system 


The VASIMR system is a high-power-density magnetoplasma rocket, which is capable of real-time 
exhaust modulation for optimum performance. In the system, shown in Figure 1 , hydrogen plasma is 
created, heated, and expelled through an open-ended magnetic configuration to provide modulated 
rocket thrust. The general magnetic structure of the VASIMR is that of an asymmetric magnetic 
mirror, comprising three linked magnetic stages. The present system preserves azimuthal symmetry. 
The magnetic configuration of the VASIMR — combined with its approach to plasma generation, 
heating, and acceleration — results in a unique engine whose theoretical performance far exceeds that 
of present-day rockets [10, 11]. 




Figure 2. Projected design of the first spaceflight low-power VASIMR thruster. 

A simpler variation of the VASIMR system is shown in Figure 2. In this configuration, the exhaust and heating 
cells are merged to provide a simpler geometry, albeit with some limitations in performance. This simpler 
VASIMR will be used in the first flight experiment envisioned to demonstrate this technology. Denoted the 
“radiation and technology demonstrator” (RTD), this experiment is currently the subject of considerable 
theoretical and experimental study. 

An experimental evaluation of the VASIMR performance is being conducted in a laboratory' device which 
uses simple liquid nitrogen- cooled copper magnets. These units are integrated into a vacuum structure by 
placing them in stainless steel enclosures, which can withstand the atmospheric pressure forces. The coils are 
attached to a structural spacer, which provides diagnostic and view port access. While this magnet assembly 
would be unsuitable for spaceflight, the resulting magnetic field closely resembles that of the actual RTD flight 
experiment and can be comfortably studied in the laboratory. Figure 3 shows a composite trimetric view of 
the present experiment configuration with an actual photograph of the device. 
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Figure 3. Trimetric view with an actual photo of the present VASIMR configuration at theASPL. 


3. The fundamental magnetostatic problem 

The magnetostatic problem is a steady state case of two vector Maxwell equations: 

Vx-B = /)! B=V*A, (1) 

where B is the magnetic induction vector, ju\s the magnetic permeability, p\s the current density and A is the 
magnetic vector potential. This paper is devoted to the fast and accurate calculation of B for a specified 
azimuthal current density p\n the special case of cylindrical symmetry and constant magnetic permeability p 
In that case, the magnetic vector potential A , written in the cylindrical coordinate system (r , 0 z), has only an 
azimuthal nonzero component: A = (0, Afr, z), 0) and the problem (1) can be rewritten in the following fonn: 

-r- — - — — = f(r,z) = jtrp, (r,z)en, (2) 
or r or oz~ 

where u(r, z) = r A fr, z) is the magnetic flux, shown in Figure 4, and /2is a computational domain. 
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Theoretically, the azimuthal component of the magnetic vector potential for the magnetic field 
generated by a coil with arbitrary cross section can be computed by numerically integrating the 
following analytical formula over the coil cross section c [12]: 




K-E 

c 

-V Z ) 

- 


dr'dz'; 


s' =■ 


4rE 

(r + r' ) 2 +(z + z' f 


( 3 ) 


where K and E are full elliptic integrals of the first and second kind, respectively [13], If the number 
of grid points and the number of numerical integration points are large, then the computational time 
can be excessive. In order to reduce computational time, equation (2) is solved in a gridded domain 
using finite difference (FDM), finite element (FEM), or finite volume (FVM) methods. In some cases, 
numerical integration of equation (3) is still useful for calculating the magnetic potential over a much 
smaller domain, such as for determining boundary conditions. 


4. Solution strategy 

In order to solve the magnetostatic problem of equation (2) in a computer we need to discretize it from 
its partial differential form into a system of algebraic equations using an approximation technique. 
Furthermore, the computational domain, shown in Figure 5, must be discretized also into a mesh of 
nodal points. Let us take on the domain discretization first. 


4.1 Domain discretization 

For simplicity, we consider a single coil with rectangular cross section c = {R; < r < R 2 , Zi < z < Z 2 } 
in the rectangular computational domain Q - {0 < r < R 3 , Zo < z < Z 3 }. Since the coil is symmetric 
with respect to z = 0, the computational domain can be reduced by half to Q = {0 < r < R 3 , 0 < z < Z 3 } 
with a boundary condition du/dz = 0 at z = 0, based on symmetry. Another boundary condition, u = 0 
at r = 0, follows from the definition of u(r, z). 
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Following the geometry of Figure 5, the above domain can be discretized as follows. Let h and h be 
mesh steps for the mesh: 


r 0 = 0 
r,=r 0 + ho 

r 2 = n + h i 

r 3 = r 2 + h 2 
r N r = r N-l +h N-l 


z 0 =0 

Z 1 ~ z 0 + fl 0 
Z 2 = Z 1 + hj 

z 3= z 2 + h z 2 

Z N, =z ^-/ + ^-7 


Finally, let us define a new array u,j such that u tJ = u(r„ zj) at each node of the rectangular grid. 



4.2 Equation discretization 

After the computational domain is discretized, the second step is to approximate equation (2) into a set 
of algebraic equations with respect to Uy. We explore two techniques to do this. The first one is a 9- 
point, finite difference, fourth-order on uniform grid (9FD4U) approximation. The complete details of 
the derivation are described in Appendix 1 . The result is the following system of algebraic equations 
with respect to uy: 


P°ij U ij ~ PhjUi-l.j ~ P 2 iJ U i,j-l ~ P 3 ij lt i + I,j ~ P 4 ij U i.j+l 
- p6 v u i _ u _, - p7 g u^,j_, - p8jjU i+ ij+i = f0 9 . 


~ P 5 y u ,-i.j+, 

i = 2,...N r - /, j = 2,...N Z -I, 


( 4 ) 


where the coefficients and right-hand side formulas are given by: 
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8 2 r 

P°i = I pk, , p2j = —j—j 2 
k 3h z ~ ( 4fj - h r ) 


( ,2 , 2 \ 

5 S— 

- K 2 ' 2 






2 2 


p3j = 


Sh/-h / , „ 1 

— = M+/> p 7 /= 


6h r 2 h z 2 r i+l/2 


1 2 r; 


i+1/2 


1 ] 


K 2 K 2 


= p8j = p5 i+I = p6 i+I , (5) 


f®ij f\j j 2 


(fij ' ( fij fi+I / 2 ) + 2 fij fi,j+l 

r i-l / 2 r i+l/2 


Figure 6 shows stencils for the 9-point and 5-point schemes described in this paper. 



A similar 9-point finite difference scheme for the Poisson equation in Cartesian coordinates is well 
known [7-9], However, such schemes have not been investigated thoroughly for partial differential 
equations in cylindrical coordinates as in our case. 

The efficient use of this method has some limitations. For example, its accuracy is greatly reduced 
when using a nonuniform mesh. Also, the grid aspect ratio (h,/hj must be close to unity as described 
in equation (19) of Appendix 1. Further still, use of a uniform mesh requires the computational 
domain size to be comparable to the characteristic dimension of the magnet. In such cases it may be 
necessary to calculate non-zero boundary conditions, using analytical formula (3). This may not be 
convenient, when computing field values over much larger domains or when computing field 
contributions due to internal plasma currents. 

The generalization of the high-order discretization using finite volume approach is described in 
Appendix 2. Appendix 3 discusses fast and accurate ways to derive boundary conditions by numerical 
integration of analytical formula (3). The results of the investigation of the 9FD4U method will be 
presented below after we discuss a second approach. 
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While the 9FD4U method has grid uniformity limitations, a simpler scheme using a 5-point stencil is 
available. This method does not suffer from previous mesh uniformity restrictions and yields faster 
solutions due to the reduced computational load. These benefits, however, come at the expense of 
accuracy. We have compared the relative advantages and disadvantages that each method offers and 
present these results in tabular form in Table 1. We now describe the 5-point scheme in detail. 

We will use the name “5FD2N” to refer to the 5-point finite difference second-order on nonuniform 
mesh technique, which can be presented by the following system of algebraic equations with respect to 

u u- 

pOjjUjj - pljjUj_]j - p2jjUj j_ } - p3jjU i+l j - p4jjUj j + ] = fOy, i = 2,...N r -I, j = 2,...N Z - 1, (6) 


where the coefficients and right-hand side formulas are given by: 


4 h'j—i + h~; 

pOjj = I pkjj , p Ijj = 77- = p3i-ij . 


n-i/2 h i-i 


K-i + K _ 


p2jj = = p4 iH , fOy = 

n h j-i 


(h[-, + h[)(hj_j+h Z j) 


( 7 ) 


2 n 


fij- 


The scheme has been previously used [2]. One advantage of this method is that the matrix of the 
system for any mesh steps is less dense than one for the 9-point scheme and the system (6-7) has a 
solution for any grid. Another advantage is that using a nonuniform mesh allows the use of a large 
computational domain without sacrificing the minimal mesh size and keeping computational error 
small. Furthermore, using a nonuniform mesh allows the use of very large domains. This in turn 
permits the boundary condition w = 0 to be used instead of the more computationally expensive 
numerical integration of formula (3). 


Both numerical approximations are compared below. 
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Features \ Methods 

9FD4U 

5FD2N | 

1) Approximation 
accuracy 

fourth-order: k = 4 

© 

second-order: k - 2 

gj 

2) Sparsity of the 
matrix 

9-diagonal 

© 

5 -diagonal 

© 

3) Mesh requirements 

Required uniform mesh with 
h/h : in some range. 
Composition of uniform 
meshes introduces 
increasing error 

© 

Nonuniform mesh works as 
good as uniform, when | h— 

hi+i\-0(h) 

© 

4) Computational 
domain 

Because of 3), domain 
cannot be very large 

© 

Because of 3), mesh can be 
large 

m 

5) Boundary conditions 

Because of 4), BC needs to 
be calculated using 
numerical integration of 
analytical formula 

© 

For large computational 
domain, BC can be zero 
without introducing big 
error 

© 

6) Iterative solver 
performance (number 
of iteration) 

Because of 3), matrix does 
not have big condition 
number, hence the solver 
needs fewer iterations 

© 

Because of 3), the matrix 
may have a big condition 
number, hence the solver 
needs more iterations 

© 

7) Number of 
arithmetic operations 
per iteration 

Because of 2), each iteration 
requires more arithmetic 
operations 

© 

Because of 2), each iteration 
requires fewer arithmetic 
operations 

© 

8) Applicability to 
currents in a complex 
geometry 

Because of 5), it is very 
difficult to use the method 
for non-rectangular coils or 
plasma currents 

© 

Because of 5), there is no 
problem using the method 
for non-rectangular coils or 
plasma currents 

© 

9) Inclusion of the 
current source in the 
domain 

If the currents source is 
inside the computational 
domain, the error may be 
large 

© 

Inclusion of the current 
source into the domain does 
not increase the error 

© 


Table 1. Comparison of the 9-point finite difference fourth-order on uniform grid approximation with 5- 
point finite difference second-order on nonuniform mesh approximation. Symbols <8> and © are used to 

show pros and cons for each scheme. 
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One of the important features of approximation methods, compared in the Table 1, is accuracy. In 
general, accuracy refers to how close the numerical solution approaches the exact one. We choose to 
measure the accuracy by bounding the difference between these two solutions to some constant times a 
power of the step size. This is expressed by 

max | Ujj - u a ( r h zj ) |< C min( h) k = 0( h k ) , (8) 

where the analytical solution if can be calculated from formula (3). The value of k refers to the order 
of the numerical error, or order of accuracy. The higher the value of k, the more accurate the solution 
becomes. In order to find a value of k, one can proceed with the full numerical and analytical 
evaluation of the left-hand side. However, it has been shown [2] that the k value can be obtained 
without requiring this lengthy approach, by using a qualitative analysis of the numerical 
approximations. Using this method, the value of k is found to be k - 4 for the 9FD4U scheme. For 
4FD2N, if a quasi-uniform grid is used, i.e. | h,- - h,.,\ = 0(h 2 ), then k = 2. k = / for both schemes, 
when an arbitrary uniform grid is used. 


4.3 Solving the algebraic system 

Both finite difference schemes represent the partial differential problem (2) by a system of linear 
algebraic equations. The 9FD4U algebraic system is presented in (4) and the 5FD2N in (6). The 
resulting linear algebraic systems can be represented by a matrix of the form: 


Ah Uh — ft- (9) 

Such algebraic systems can be solved either by direct or iterative methods [2, 14], 

The resulting linear algebraic systems, presented in this report, are solved by the fast iterative implicit 
incomplete factorization method [14 - 15]. The stop criteria for the iterative process of calculating u/, 
is the following inequality: 


fh A h U h 
Ifh - A h u h ° 


< £, = 10 




( 10 ) 


The number of iterations n(£j depends on the computational domain and grid and normally is between 
20 and 40. 

5. Numerical experiments for a single coil 


5.1 ID tests 

A number of numerical experiments were performed to test the accuracy of the numerical methods. 
The first and simplest corresponds to the infinite solenoid representing the ID solution. The solution 
results are shown in Appendix 4. The error for the 9FD4U scheme is virtually zero. For the 5FD2N, 
the error is very small and decreases as h 2 . The second and more interesting set of tests is shown 
below and corresponds to a single coil of finite size. 
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5.2 Coil complexity 

The complexity of solving a magnetostatic problem for a finite rectangular cross-section coil is affected 
by the coil aspect. This is so because the mesh must be compatible with the coil shape. From the 
mesh-generation point of view, it is difficult to discretize the computational domain for a very thin, 
very narrow coil or a coil with a very small hole, as illustrated in Figure 7. In general, a nonuniform 
mesh is preferable when dealing with a difficult coil. 


Difficult coil 


Easy coil 




Too narrow 



Figure 7. Difficult vs. easy coils from the mesh-generation point of view. 


5.3 Deriving an exact solution and calculating the numerical error 

The first test was considered for the “easy” coil “cl” with a cross section c = [0.5, 1] x [-0.5, 0.5] and a 
current density corresponding pp = 1. The solution for this problem is plotted in Figure 8. 



Figure 8. A magnetic field near an “easy” coil. Left: magnetic field lines (contour lines for the magnetic 
flux u(r, z)). Right: axial magnetic field, calculated numerically and analytically along r = 0. 
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To study the accuracy of the two methods, the numerical solution was compared with a known 
analytical solution calculated on the symmetry line r = 0 by the formula 


B° z (0,z) = 


(z — Z,)ln 


R 2 + yj ~R~ 2 + (z — Z ; ~) 2 
R I +y[R] ~+(z-Zr)* 


— (z — Z 2 )ln 


R 2 + yj~R ] Hz — Z 2 )‘ 
R.+^Rj +(z-Z 2 ) 2 


( 11 ) 


where current density p is assumed constant in the coil. The following maximal and point-wise 
relative errors on the symmetry line were observed: 


<5 = ma x5 b (z); d B ( z j) : 




max 

j 


B"(0.*j> \ 


( 12 ) 


where B h (0,z) is the numerical solution calculated from the numerical magnetic flux solution u at the 
symmetry line r = 0, according to the scheme described in Appendix 5. 



Figure 9. Magnetic field B, (red line) and relative errors Sr of the magnetic field Bfor the single coil “c 1 ”, 
calculated on the symmetry line r — 0 by 5FD2N (blue lines) and 9FD4U (green lines) methods using three 


different mesh sizes: 64 x 64, 128 x 128 and 256 x 256. 
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5.4 Finite difference grid 

The first comparison of accuracy of the two methods considers the dependence of the relative error on 
the grid size. This is shown in Figure 9. 

Two families of error plots corresponding to 5FD2N and 9FD4U methods are shown in blue and green 
color, respectively. In each family three mesh structures were used, including 64 x 64, 128 x 128 and 
256 x256 grid points. Several features of Figure 9 merit discussion. 

The second-order non uniform-mesh method 5FD2N is generally less accurate but has more constant 
error along the domain. The fourth-order uniform-mesh method 9FD4U produces more accurate 
solution but shows diverged error over the computational domain. 

In all cases, the error decreases with increasing grid size in accordance with theoretical predictions. 
For example, the second-order 5FD2N method makes the maximal relative error <5 decrease by factors 
of 4 as 0.003, 0.0008 and 0.0002, when the grid step size h is divided by 2. On the other hand, the 
fourth-order 9FD4U method makes the relative error <5 go down by factors of 16 as & 10' , 5. 10 , and 
3.-10' 7 , when the grid step size h is divided by 2. The strong discontinuities in the error plots are due 
to several factors, such as the behavior of the analytical solution and its derivatives, the choice of grid 
and the method used. The use of a fully adaptive grid will tend to smooth out these irregularities. This 
is beyond the scope of this paper. 

One can see that, to meet the current experimental requirements on accuracy of 10' 3 , the scheme 
5FD2N should be used with a 128 x 128 mesh. For the same requirements, the 9FD4U method can be 
used with a 64 x 64 mesh. This shows an advantage of the 9FD4U (uniform) method when the 
magnetostatic problem is solved for an easier coil. 

Mesh configurations for both methods, as applied to an easier coil, are shown in Figure 10. The 2D 
nonuniform mesh in Figure 10 (b) is a product of two nonuniform ID meshes. Each of them is fine 
next to points of singularity: edges of the coil and the symmetry line r = 0. It was created 
independently in the r and z directions, as described in Appendix 4 for the ID case. 



Figure 10. An example for the uniform (a) and nonuniform (b) finite difference meshes for “cl" coil. 
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5.5 Domain size and boundary conditions 

In the above calculation, the size of the computational domain: Q = {[ 0 , 8 ] x[- 8 , 8 ]} was chosen to be 
about ten times the size of the coil. The requirement is particular to the magnetoplasma thruster 
design. If the magnetic field needs to be calculated only near the coil, the computational domain can 
be much smaller, simplify ing the problem considerably. 

Choosing an appropriate domain size is related to choosing appropriate boundary conditions. As 
mentioned earlier, there are two ways of handling boundary' conditions. First, a semi-analytical 
formula given by equation (3) can be used to compute the field value at the boundary for both 
methods. Second, the magnetic potential is assumed to be zero at the boundary, which is a good 
assumption for a large domain. 

In the first approach, formula (3) requires the use of a numerical integration technique described in 
Appendix 3. Convergence of the boundary condition calculation method is presented in Appendix 6. 
This method works well for both 9FD4U and 5FD2N schemes, when the current source is given as a 
coil with small cross section. 

For the large enough domain, the 5FD2N method can work efficiently with zero boundary conditions. 
This property allows it to work efficiently for problems with complex current source terms, such as 
plasma currents, when the first approach based on formula (3) cannot be used. The accuracy of the 
5FD2N scheme with zero boundary conditions is demonstrated below. The 9FD4U scheme is taken 
out of comparison there, because a large computational domain with a difficult coil requires too large a 
uniform mesh. 

Since an exact solution may not be small on the boundary points close to the coil, the use of a small 
domain with zero boundary conditions could be a source of the error. On the other hand, a large 
domain could be an error source as well, if the grid is too coarse. The optimal domain size depends on 
the coil geometry and can be found from numerical experiments. The dependence of the numerical 
error from the domain size R 3 and coil geometry for the fixed number of grid points is presented in 
Figure 11. 
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Figure 11. Maximal relative error 5 of the 5-point scheme 5FD2N, as a function of the domain size Rj. 2D 
calculations of the magnetostatic problem for “cl ” and “ c2 ” coils were performed, using semi-analytical 
and zero boundary conditions. For each Rj, the nonuniform mesh 128 by 128 was built to minimize the 

relative error for ID problem. 


Let’s analyze the error for the “easy” coil “cl” and the long-shaped “difficult” coil “c2” with cross 
section c = [0.049, 0.053]x[-0.2, 0.2], shown in Figure 11. The use of the semi-analytical boundary 
conditions provides less error, especially for the small domain size R 3 . Use of zero boundary 
conditions produces a big error for the small domain size. Also, for sufficiently large domain size, the 
error would grow, because the grid becomes coarse. The zero boundary conditions yield a minimum 
error corresponding to some “optimal” value for the domain size R3. 

Since the coil “c2” has a smaller and thinner cross section than “cl”, the optimal domain size is 
different. While the total number of mesh points is fixed at 128 by 128, the plots have minimums at R 3 
= 50, 5 = 0.0003 for coil “cl” and R 3 = 8 , 8 = 0.0002 for the coil “c2”. These numerical tests 
demonstrate the robustness and good performance of the 5FD2N scheme for solving the magnetostatic 
problem for a difficult coil. 
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6. Numerical experiments for a multi-coil system 

Our ultimate test demonstrated in this section corresponds to the multiple-coil system involved in the 
magnetoplasma thruster. This system consists of rectangular cross-section coils of different shape. This is 
why our magnetostatic solver needs to work well for a variety of coil geometries. 

In practice, the usual problem consists of computing the magnetic field of the set of various coils of different 
sizes and placed at varying distances from one another. Fortunately , the principle of superposition applies and 
the total field of the coil system is the sum of the fields of separate coils. In doing so, the size of the 
computational domain and the number of mesh-points can be decreased dramatically. As shown in Figure 12, 
a nonuniform mesh has singular points at the coil edges; hence the number of singular mesh points is bigger for 
the multi-coil system. By breaking up the magnet into single coils, the number of singular mesh points is 
greatly reduced. 



Figure 1 2. Illustration of the principle of superposition. The magnetic field for the multi-coil system 
(top) can be calculated as a superposition of the magnetic fields of each individual coil (bottom). Sample 

nonuniform meshes are shown for the 5FD2N method. 


The full magnet field of the low- power VASIMR thruster is shown in Figure 13. Table 2 specifies the various 
magnet parameters. 
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Figure 13. Magnetic field for the 14-coil magnetoplasma rocket configuration. 
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N 

R, 

Ri 

z, 

Z 2 

J 

Section 

Sample nonuniform mesh 

1 

0.1068 

0.1108 

0.02 

0.04 

2000.0 

Helicon 

magnets 

i 

t* 

i 

i 

1 

Hutu 

2 

0.1028 

0.1068 

0.02 

0.06 

4000.0 

3 

0.1004 

0.1028 

0.02 

0.12 

6000.0 

1 "* tiJi ifciiUimii ii SmJt 

4 

0.098 

0.1004 

0.02 

0.16 

8000.0 

IfftiiSflfl 

5 

0.095 

0.098 

0.02 

0.18 

12000.0 


6 

0.095 

0.145 

0.3075 

0.3525 

110000.0 

Choke 

magnet 

11 !.;■): :: iil ll illlir ::x 

H ■" •; -~ij» 

;| || 

7 

0.095 

0.101 

0.39 

0.53 

42000.0 

ICRH 

magnets 

i 

1 


8 

0.101 

0.103 

0.45 

0.53 

8000.0 

9 

0.103 

0.105 

0.47 

0.53 

6000.0 

I ■ ”"-»*■ ' r '" T " 

10 

0.105 

0.107 

0.49 

0.53 

4000.0 

11 

0.107 

0.109 

0.51 

0.53 

2000.0 



12 

0.113 

0.117 

0.53 

0.61 

16000.0 

Magnetic 

nozzle 


^ISlfiiyspa 

13 

0.117 

0.125 

0.57 

0.61 

16000.0 



14 

0.125 

0.145 

0.59 

0.61 

20000.0 




: -Sli-i ti t : -ini-ifca-;. jjArtaSSifc 


Table 2. Description of the multi-coil magnetic configuration considered for the RTD VASIMR thruster 
shown in Figure 13. Here J is the total current in each coil. The current density function p in each coil is 

computed as p = J / (R 2 - Ri)(Zi - Z,). 


The uniform and nonuniform meshes used are demonstrated in Figure 14. For the 9FD4U method, the 
computational domain was chosen as the union of two rectangles: 0 < r < 0.075, - 0.1 < z< 0.7 and 0 
< r < 1, 0.7 < z < 3. Such domain includes the area where plasma flow is observed. Since “difficult” 
coils require too-large uniform meshes, the uniform-mesh domain excludes the coils. A piece-wise 
uniform grid with 16 x 72 points is used for the first subdomain and the grid with 240 x 64 points is 
used for the second one. The 5FD5N method used nonuniform mesh with 100 x 300 grid points on the 
domain [0 5] x [-5, 5 ]. 
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Figure 14. a) Sample uniform meshes used for the 9FD4U method (top half of the picture); area around 
coils is not covered by the 9FD4U method; b) sample nonuniform mesh used for the 5FD2N method (bottom 

half of the picture). 


The accuracy of both methods is demonstrated in Figure 15. It has been found that there is about the 
same level of the numerical error for both methods, which is below the current experiment 
requirements of H)' 3 . The numerical error for the 5FD2N method with nonuniform mesh is closer to a 
constant than the error for the 9FD4U method with uniform mesh. The uniform mesh error is higher at 
areas corresponding to a larger magnetic field. Also, the uniform mesh error has a discontinuity at the 
interface boundary of sub-meshes. 
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Figure 15. Relative error of the magnetic field for the set of 14 coils calculated on the symmetry line by 

5FD2N and 9FD4U methods. 


7. Computational complexity 

The computational complexity involved in the solution of the problem generally implies requirements 
on computer memory and speed. For some applications these resources may be limited. Hence the 
proper choice of computational method is determined in the context of the full engineering problem. 
To illustrate this point, we examine the tradeoffs in complexity as functions of the computational error 
for the simple coil shown in Figures 8 and 10. 

First, let us discuss computer memory requirements. One of the most efficient available iterative 
solvers, the iterative implicit incomplete factorization method (IMIF9, [14]), requires the use of a 
memory volume P « T AM N r N-, where T AM is the number of 2D arrays used by the approximation 
method. The size of arrays is equal to the number of mesh points of N r N : . The 9FD4U method uses 
T am = 12 arrays, and the 5FD2N method uses T AM = 9 arrays. Therefore, using double precision, a 
conventional personal computer with 128 MB of RAM should be able to handle the problem as long as 
the mesh does not exceed 1000 x 1000 points. 
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The requirements on computer speed are a bit more complicated to evaluate, as they depend on the 
number of iterations in addition to the mesh size. An approximate formula for the number of 
arithmetic operations Q, required by a numerical method to get a solution with given relative error 8, is 
given by Q~ U at n(£,) N r (5) N z (8). Uat is a number of arithmetic operations used by a method, per one 
iteration, per one mesh point. 9FD4U method implies Uat = 45 and 5FD2N method implies Uat = 21 . 
n(£j) is the number of iterations required to get numerical solution with an error less than The 
iteration error £, has to be much less than the relative numerical solution error 8. One can use the 
following simple relation: e, = 8 2 . As numerical experiments show, the number of required iterations 
n(£i) is between 20 and 40 for the iteration error £, = 10' 6 . 

The mesh size for the fourth-order 9FD4U method is bounded by a numerical error as N r < 0(8 14 ), N z 
< 0(8 14 ). For the second-order 5FD2N method, the mesh size is defined by the relation N r < 0(8 1/2 ), 
N z < 0(8 I/2 ). All this makes the Q(8) dependence to be quite complicated. Rather than attempting to 
evaluate this expression, a set of numerical experiments has been performed for a simple coil (cl). 
The results are plotted on Figure 16. The number of arithmetic operations is measured in million float- 
point operations (MFlops). 
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Figure 16. Numerical error 8 vs. number of arithmetic operation Q needed 
for magnetic field calculation in cl -domain. 
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From this figure, if one is interested in high accuracy, the 9FD4U method is faster, whereas if accuracy 
is less important, the 5FD2N method is more desirable. As pointed out earlier, these results are 
dependent on the complexity of the coil structure. For a relative error close to the current experiment 
requirement of <5 = 0 . 001 , the 9FD4U and 5FD2N require about the same number of arithmetic 
operations for the “easy” cl -coil problem considered here. 

However, the second-order scheme using nonuniform mesh is a better tool for calculating a magnetic 
field for more complex problems. One such complex problem is a multi-coil system including 
“difficult” coils, when the magnetic field has to be calculated next to the coils. Another complex 
magnetostatic problem is solving for the magnetic field in the presence of internal plasma currents, 
where no semi-analytical solution for the calculation of boundary values is available. 

The total computational cost of calculating magnetic field must include the cost of computing the 
boundary conditions, when they are calculated by the numerical integration of formula (3). For simple 
coil geometries, the numerical calculation of the boundary magnetic field requires many fewer 
computations and gives much less error than finite-difference solvers. This is why the computational 
complexity of boundary condition calculation does not affect the overall complexity of the magnetic 
field calculation. 

8. Conclusion 

This paper demonstrates advanced computational techniques for an accurate computation of the 
magnetic field of an axisymmetric electromagnet. These techniques involve the adaptive and uniform 
grids, second- and fourth-order approximations, and use of a semi-analytical calculation of the 
boundary conditions. The comparison analysis has been conducted for two approximation methods: 
the fourth-order 9-point finite volume scheme using uniform mesh (9FD4U) and 5-point finite 
difference scheme on adaptive nonuniform mesh (5FD2N). The approximation method of choice 
depends on the complexity of the current sources (coils), specification of the computational domain for 
finding the magnetic field, computational accuracy requirements, and available computational 
resources. 

For achieving the accuracy imposed by our experimental requirements, both 9FD4U and 5FD2N 
methods have about the same computational complexity. If higher accuracy is needed, the fourth-order 
scheme 9FD4U will have dramatic advantage over the second-order method 5FD2N. The 5FD2N 
scheme is the best tool when the computational domain has to be of the large size or there is no semi- 
analitical method for boundary values calculation. Such a case exists when the presence of internal 
plasma currents affects the overall magnetic field. The numerical methods described in this paper have 
been applied for solving an electrostatic problem, which is a part of the extended mathematical model 
in the magnetoplasma rocket simulation [2]. 

Numerical predictions of the field generated by an axisymmetric electromagnet are also being applied 
to experimental investigations of the VASIMR system in the laboratory'. The magnetic field contours 
shown in Figure 17 are outputs of this predictive tool. Actual photographs of the plasma show plume 
divergence as predicted by the code. 
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Figure 1 7. Magnetic field in the present VASJMR configuration calculated by the numerical code 
together with superimposed photographs of plasma source and plasma exhaust. Magnetic field lines are 
plotted for the area with a plasma flow, as well as the magnetic field strength contour lines. The bottom 
part of the figure demonstrates a ID plot of the axial magnetic field along the symmetry line. 


The 5FD2N method is currently used operationally to design an experimental configuration of the VX- 1 0 
device at the ASPL. For given electromagnet geometries and current values, the magnetic field is calculated 
and plotted (Figure 17), as well as contour lines for the magnetic field strength. The present code allows the 
determination of the electromagnet currents needed to generate a required magnetic field. Another important 
application of the magnetostatic solver based on the 5FD2N method is a self-consistent simulation of plasma 
flow in the VASIMR. 
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Figure 18: Local numbering used in the compact 9-point finite difference scheme (17) and two boxes 
around mesh node used in finite volume discretization (21). 

Appendix 1: Compact finite difference scheme of the fourth order 

Our goal is to derive the compact finite difference approximation of equation (2) of the fourth order 
at the uniform grid with mesh steps ft- = h r . h z - = h z (here “compact” means that only the closest to 
the central point mesh points are included into the stencil for approximation). To get this, we define 
the scaled one-dimensional (ID) differential operator 

t 1 r dldu 

L t U = -L r U — - — - — ( I'D 

r or r or 

and its three-point approximation with the following property: 

dSkj = ^( = (L r u - %.L r ri r u)ij + 0(h 4 ), (14) 

l} r \ D— 1/2 r t'+l/2 / 1Z 

where r i±1 ^ 2 = (r,- + r t ±i )/2. The validity of this relation can be checked by Taylor expansion. 

From (14) and similar definition of another ID difference operator 

( Lz u )i,j - J~2 ( 2u *,i - “ w f,j+i) = (L z u - Llu)ij + 0(h 4 ) (15) 

w r e have the following nine-point, scheme of order 0(h 4 ): 


LU)u = rl h T u + L n z u - —Ah 2 r L h r L h z u + h'iL n z L n r u 


.2 T h T h. 


= fij - —(rh 2 r L h T .f + hlL*f)ij + 0(h 4 ). (16) 

The validity of this relation is proved by the direct substitution of the expansions (14), (15) in (16) 
and using the original equation (2) for the terms of the second order. The function / is required to 
have bounded fourth derivatives. 

Let us introduce the local numbers for the nodes of the nine-point, mesh stencil in accordance with 
Figure 18. Then we can write the resulting grid equation of the mesh point (?. j) in the form 


pOijiio ~ = JOi 
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where the coefficients and right-hand side are given for r, : > 0 by the following formulas derived by 
multiplication of ID finite-difference operators: 


P^i,j 

P3«j 

P2«J 

P"i,j 

/Oy 


fc=l 

1 


h 2 

(5 1 2 1 


1/2 l)2 z 

2 p /g ^ o h r , _ . 

3ft*(4r?-ft*) (5 h 2 3 2/f P ,J 

77 + 77) = P%i,j — 1 = PX+1J— 1 = pbj'+ij. 


12r t -n/2 h 2 h 2 
fi 


,3 ~ tx[ — ' — (/•',; - fi-1,3) d (/*,i _ /i+ij) + 2 /tj /<j+i]- 

12 r t _!y2 


'i+1/2 


All the coefficients (18) are positive under the conditions 


2/7 < h 2 //* 2 < 5. 


( 18 ) 


(19) 


If we have the Dirichlet condition on all of the boundary OQ of the computational domain u\ du — g(r, z) 
(the simplest case is g = 0), then an algebraic system of linear equations (17), (18) can be rewritten 
in vector-matrix form 

A h u h = A, Mfc = fh = {Aj}- ( 2 °) 

where is a symmetric and positive definite matrix for an arbitrary ratio of the niesli steps. Ah is 
a Stieltjes type matrix, i.e. it is monotone (A^ 1 > 0, see [14] for example) only under the conditions 
(19). The convergence rate O(h^) of the error of this scheme can be proved in the Euclidean or 
maximum norms by usual techniques, like it was done in [7, 8]. 

The generalization of the high-order discretization using finite volume approach is described in 

Appendix 2. 


Appendix 2: Finite volume approximations 


Using direct finite-difference approximations, it is impossible theoretically to get high-order accuracy 
if the functions u and / don't, have smooth enough derivatives. 

So, our second approach consists of constructing the finite volume scheme for the piece-wise smooth 
function /. The approach is based on the approximation of the integral equation obtained by formal 
integration of the differential equation (2) and coincides with the fourth-order equation (16) in the 
case of a uniform grid and a smooth function /. 

A similar finite-volume scheme for the Poisson equation in Cartesian coordinates can be found in 
[8, 9]. Such schemes are not investigated thoroughly for partial differential equations in cylindrical 
coordinates. 

In accordance with [9], we introduce the small and big boxes around the (i,j) mesh-point: 


Vi,j = (r,- - /ij_i/2 < r < + fii/2, Zj - fij-i/2 < z < zj + hj/2^j, (21) 

Vi,j = (d-1 < r < r t+l . Zj- i < £ < 2j+i} 


and denote their boundaries by S t j and S itj respectively, as shown in Figure 18. 
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Let us integrate by parts the left side of equation (2) over the volume 



u + L z u)dY - 


Zj + l/2 
/ 




Z J~ 1/2 


r t41/2 



T i- 1/2 


1 Ou 
r dz 



and consider the linear combination with a weight parameter 0 < ir < 1 of the integral relations for 
the small and big boxes 

lij = u- j r n dS + (1 - tr) J I n dS - w j ~drd~ + (1 - «') J '-drdz. (22) 

s i,i S u Vij v tiJ 


1 () u d it 

Here I n is the projection of the vector I = (/ r , I z ) = — We assume also that some 

jumps of function / are possible only at the grid lines r = r*, z = Zj, 

To approximate equation (22), we introduce sub-volumes V-j, I = 1,2, 3, 4 which are the 

result of subdivision of Yij and V{j by the grid lines into four parts: 


4 



(=1 


4 



?= 1 


We define corresponding fluxes ijj through the sub-volume surfaces Sjj and 



4 



1=1 


(23) 


and write with the help of (22) the necessary relation for / — 3 in more detail (right upper sub-volumes 
for r > r*, z > zj): 



Additional terms with the integrals over the “interna]” boundaries between sub-volumes are included 
here to improve the approximation of the derivatives. It is important that in the sum (23), additional 
terms annul because of the continuity property of the solution and its fluxes. 

Then by applying the simplest quadrature formulas and the linear interpolation of the terms under 
the integrals, we can write the flux l\^ as some linear form of four values u\: 

= 0 + p03ij«3 + p0ii,jU4 + p08jjW8- (25) 


If the fourth-order vector uf 3 ) = (w 0 , w 3> M 4: w 8.)? whose components are the values of the solution in 
the vertices of the finite volume I,j\ and vector / ^ = (I^\ l[ 3 \ 1^), whose entries are the 
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“partial" fluxes around corresponding nodes (for example, if^ - ff9 from (24)), are introduced, then 
one can write in vector form the following relation between lliem: 

/pOOij p03ij pO-iij pOSij 
7 ( 3 ) _ ,( 3 )-( 3 ) 4 ( 3 ) _ P30ij p33ij p31,j -pMij 

piOij p43ij p44,j p4Hij 
V pSOij pS3ij pS p 88 tj 

Here the fourth-order symmetric matrix (pfc/.j = p/Aqj ) can be called the local balance matrix 
by analogy with the local stiffness matrix of finite element methods, introduced in [16]. 

It is easy to show that the elements pk hJ of the “global" balance matrix A h from system (17), (20) 
can be written via the entries of local balance matrices by the formulas 



p0,-j = pOOij + p33,_i j + p44i t j-i + p 88 ,_ij_i, p3 tJ - — pUj — p03ij + p08j (J _i, 
p4-j = p2i tj = pOlij + p38j_i,j, pSij = p7ij = p5ij = p6,j = p08 i>} . 


In practice, instead of implementing these formulas, it is preferable to assemble the global matrix 
Ah in an element-by-element technique which can be presented by the following pseudo-code for each 


sub- volume \)j: 


pO t j := p0t,j + pOOjj, p0j'+i,j pOj+ij 4- p33,j, 

p0,j + i := pO tJ+ i + p44,j, pOj+ij+i : = pO*i,*i + P 88 «,.r 

p3,,j p3i,j + p03,j, p3i,j+i := p3,j+i + p48^ v (28) 

pUj := p4 t> j + p04,j, pl.+ij := p4<+i,i + p38,j, 

P 8 ,j := P 8 i,j + p08qj, p7,-,j+i := p7»',*i + p34;,j. 


It is possible to use different approximations of the integrals in (24). We apply quadratures which 
follow from the linear interpolation of the functions under the integral over the surface of l' t J : 


— if j — w s , [3( j w i+i j ) A Ut,j+i u *i,*i] 

’ J l 8 /! i r i+l /2 

+ 15 $ [(^) ( “ itu _ “ i+w+l) + ( 77i+,/2 " 7 ^) ( “ iJ " “” i+1 ’] } m 

+(1 - w) 77— 2 («»,/ - «* 1 ,j + «ij + 1 - Mt+lj+l) + _ Ui ’j +1 + ? '*l,i - «*l,J+l) • 

.^j r *l/2 j 

Here the values 7 i+1 , 2 are defined from the condition that the entries p/qj of the global matrix Ah- 
which are assembled via the elements pkh,j of the local balance matrices by the formulas (25) - (28), 
coincide with the coefficients (18) of the compact finite difference scheme of order 0(h 4 ) in the case 
of uniform grid and w = 16/15. This condition has the recurrent form 

h] + /({_! ^N r +i ■ y I (3Q) 

7i-l/2= 7*1/2, /JVr+1/2- » l ' iU > 


lr i + n i-l . . _ "M+l ; _ V 

7,-172 = —-7*1/2, /AW 1/2 - rjVr+1/2 * - 

and provides the following formulas for the entries of the local balance matrix: 

1 [7 5 w\ h z j w /([] f 15uA 7, 

p03 i , j = p\8i,j - — — (l“ 8 ) ft r 16 IP V 16 J 


1 _ 

8 / ftj 1 6hj 


(-04,, = = (l - ^) 7f 


pORqj - p34,j — 


, _ Ie) 2 + +(, 

s j ir> j \ 


15u>\ 

7*1/2 

16 ) 


7 iv 

A hz A 

~ ~8 

J hr \ ’ 

15w\ 

7*1/2 

16 J 

^2 


pOOij = p03,j + pOljj -f p08qj — p33( t j — p44 tJ - pS 8,j. pklij — plki j. 
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The approximation of the right-hand side of equation (22) can be made in a similar cell-by-cell 
approach, so the value / 0,j from (20) is assembled as 


fn . . /-( 3 ) I n j ) i f( J ) , i fW 

1 ' ./8,4— lj — 1' 


r(3) 


(3) 


f(3) 


(31 ; 


(3). 


Here the following ‘‘local" terms are defined from the numerical integration on the sub- volume lb : 


r (3) _ r M _ « ~ b w jimj - Jj+iPi± i,j 

i “ [v 1 ^ W,, 

,(3) [, , 3u’ 8 - 7t» , r.+iPi+i - 

= [(1 - — ) / > «+ 1 

= [ (1 “ T )/),J+1 “ 16 r <+1/2 

,(3) ( M 3tr N 8 - 7w^ r-, + ip, + ij + x - r t 7>,j + i t ^ 

~ [(1 ~)/ > «+lJ+l ( r-.. ^ Pi+1,3+1 


4 - ig~ ( 


)]*<*?. 


3tP^ 8-7tt7i +1 p w -Wj , „ . ^ Ur^ 

-r)/>i+ij 77 ( + /h+ij “ A+i,j+i ) pV ? j 

4 16 r i+i/2 

3u\ 8 — Iw t Ti+ipij+i — r*+i/>t+i j.fi l ^ 

77 ( r Pij + l 


*+1/2 


(32) 


(34) 

(35) 


where pij = /(r,-, 2:j)/r,- at the line of jump corresponds to “own” sub-volume. It is easy to check 
that formulas (31)-(35) provide for u? = || the term /O^yfrom (18) in the case of a uniform grid. 
One remark: in fact coefficients in (18) in Appendix 1 differ by the multiplier (1 — |u?) (h? + h ) 
(h z - + hj_ x )/4, which is the “weighted” cell volume. 

These formulas provide the symmetry of the local balance matrix and global matrix A h,. The 
values *fi±i /2 are positive for the quasi-uniform grid. For the general nonuniform grid, the approx- 
imation of this scheme is of order 0(h) only and the matrix Ah is monotone under strong enough 
conditions for the ratio of the mesh steps (19). 


Appendix 3: Computation of boundary conditions 


It is known [12] that the magnetic field of a single current ring for the unit magnetic permeability is 
given analytically by the following formulas for the axisymmetric azimuthal component of the vector 
potential A and r- and ^-components of the induction vector B : 


A e (r,z) = 
B r {r, z) 


4 J 


T 

B Z {T,3) = 



(r+r') 2 + (~ — ~') 2 ' 


2J 

V+'+ r)2 + (;-?)* 

2J 

\f(r' + r) 2 + {z - z 1 )^ 


(V) 2 + r 2 + (z - z') 2 
( r'-r f + (z - z') 2 

(r') 2 - r 2 - (z - z') 2 
(? - r) 2 + (z - z') 2 


E - Aj , 
E + h] . 


(36) 


(37) 


(38) 


Here r f is the ring radius, J is the ring current, z* is the coordinate of its center in cylindrical coordinates 
r, z and A’, E are full elliptic integrals of the first and second kind respectively (see [13]), for which 
the following series are valid: 


;*'-5 

— E = -K 

7T 7T 


1+( 2 


« 2 + 


1-3 


s 4 + 


Ks* + 


2-4 

2 / 1 - 3 N 2 


1- 3-5S* 

2- 4-6 


, 6 + ... 


3 \2 • 4 


2 3 

" + 5 


\2 • 4 • 6/ ,\2-4-6-K) 
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By the definitions of .s, E and A , it is possible to show that .1(0, z) — B r {0,z) — 0. 

The computation of the magnetic field generated by a coil with arbitrary cross section can be made 
by numerical integration of analytical formulas (37), (38) over the coil cross section 

c={0<Ri<r'< R 2 , Z x < z 9 < Z 2 }. (39) 


The following formula, which can be obtained from (36), is used to compute the boundary values 
of the solution: 


u{r,z) 


“2 1 S 2 2 

pjf \/(r + r') 2 + (z - z') 2 M(s 2 )dr , dz', M(. s 2 ) = -(1 - - —E, 

Ri Zi 


(40) 


where p is a uniform current density and the coil cross section is the region R i < r' < R 2 , Z\ < z' < Z 2 . 
For the underintegral term M (s 2 ), the following series can be derived via the expansion for full elliptic 
integrals K, E: 

OO 

= E<v ,4+2p - ( 41 ) 

P = 0 

Instead of M{s 2 ) we use a calculation of its approximation for s < 1: 


Mis 2 ) = s 4 


1\ 2 1 
2/ 8 + 


1 • 3 
F~4 


5 + 

12 


i^y^ 

2-4-6/ 16 


M n ( i 


= E< 

P = 0 


t .4+2 p 


(42) 


To understand the convergence rate of M n (s 2 ) when n — ► oc, the values of 


= 1 + 


{r-r’)* + {z 
4 rr f 


J\2 


-l 


can be estimated for (r, *) at the boundary of the computational domain (r — /? 3 or 2 — Z 2 ) and 
(V, z r ) from the coil cross section (39) by the following inequality: 


s 2 < max 



(i?3-i?2) 2 \ 1 

4 R 3 R 2 ) 



(_Z3-_Z lf\ M 
AR 3 Ri ) J 


For example, if R 3 > 3f? 2 and Z 3 - Zi > 2R 2 , then s 2 < 0.75. So, if the domain boundary is away 
from the coil, then the value of s 2 decreases and the error of approximation of M(s 2 ) by M n (s 2 ) for 
the fixed n is improved. More exactly, the cutoff error e n can be estimated by the following relation: 


— U 


ji/„ 


£ «, 

p = 7l + l 


F +2p < c„ +1 .s 4 + 2 (” +1 )^s 2 ' < c„ +1 .s 4+2 (” +1 )/(l -s 2 ). 

1=0 


Our computational strategy is a mixed one and consists of numerical integration with respect to 
r' by the simplest central rectangle quadrature rule 


/ -R2 m 

*J n { v')dr' + C~n — pile E u; u( r 'l ) + 

/— i 

Ri 


+ c! 


(i) 


r,' = «!+(/- j )/‘c 


K = 


Ri — R\ 


m 


(43) 
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or by the Simpson’s quadrature 


u(r, z) = PjMr?) + -U’„( r") + 2^ n (r") + . . . + u,„(r'' m+1 )) + c n + rf = i?i + (4-1) 


or by the Gauss quadrature 

u( 7 q z) — ph c + • • • + fm^n(^m)) T T 


(45) 


where 


^2 

-’„(>') = j y/{r + r') 2 + (z-z')Hl n {s 2 )dz', 


z 1 

/?2 Z2 


e n = p J J \J (r + r ') 2 + (2 — Hdj*t n dr'dz' < p\f(Jh + 7? 2 ) 2 + (Z 3 - Zi) 2 (/? 2 - i?i)(Z 2 - Zi)e n 

Zj 

with the given integer m. Here i = 1, 2, 3 are the errors of different quadratures. The total error 
of such approximation is 

e = e'n + = 0{h1 + S 2n+6 ), i = 1,2, 7= max {s(r, r\ z, 

(r . « ) € C, 

(r, z) € 

where 7 = 2 for quadrature (43) and 7 = 4 for (44) and can be made as small as necessary with 
sufficiently large integers m and n. The optimal Gauss quadrature points r*. and the coefficients Ck 
in (45) are defined in [17], for example. The error of the Gauss quadrature is estimated by the 
inequality 

m ~ [(2ro)!] 3 (2m + 1 ) Ri<r’<R 2 \ {dr f ) 2m 

As for the error c n , for s 2 < 0.75, for example, we have C 23 < 4.2 • 10 -9 , C 33 < 9.5 • 10” 12 . 

To compute ^ n (v r ) we use the auxiliary recurrent relations for the following integrals: 




<t>2p+3(z, t) 


dx 

(Vx 2 + /)p’ 

1 ( X 

t{2p+ 1) l (v^2T7)2 P +1 


dx 


( y/x 2 + /) 3 / \/.r 2 + / ’ 

~t~ 2p<^ 2 p+l) , P — 1)2,... 


(46) 


If we denote now 


<f> p (r, r\ z, Z lr Z 2 ) = <f> p (z - Z 2 , (r + r') 2 ) - <? P 0 - Z,, (r + r') 2 ), 
then it is possil)le to write the resulting formula to compute u; n (r[): 

u<’„(r{) = 57 Cp( 4 r / 'r) p+ 2 0 2 p+ 3 ( r ’ r (- -• z i- Z 2 ), (47) 

p—0 

where coefficients c p are defined in (41). 

The stable computation of «; n (r[) is held by using the recurrent procedure. Namely, we define the 
new functions 

d>2 P +3 = (4r r') p+ 3 / 2 5 2p+3 . p = 0, 1 , . . . , 
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which arc computed from the relations 


7 1 
<t> 2 p +3 - ^ 2 p + l) ( j 



ft + 


JEM 

4 TT f 


2p+l 


V-irr’ + 2])6 2 p+i 


) 


( r + r') 2 

Jr? 


and calculate «; n (rf) by the formula 


M?i) 


^n-l( r \) + ^2n+3 


— IT 

8 + 4(« - 1 ) iix 


/ 2m - 1 

\ 2?n 


2 


The presented formulas provide high accuracy of computing the boundary conditions for increasing 
integers m, n if the value of s 2 is not very close to unity. 


Appendix 4: ID test 

This section describes convergence of finite difference methods of the magnetostatic problem in a 
case when the solution depends only on one (radial) variable. Then we can present an algorithm of 
generation of ID nonuniform meshes. 

Consider a coil in the shape of the infinite solenoid with internal and external radii R\ and R 2 
correspondingly. Assume that the coil has an azimuthal electric current with unit density, and the coil 
is placed inside an infinite conducting cylinder with radius 7? 3 . The magnetic field for such system is 
one dimensional, B = (0,0, B z (r)) and has the following analytical expression: 

f R 2 — R\ - 2d, 0 <r< f?i, 

B z (r) = J? 2 - 2 d - r, R\ <r< f? 2 , (48) 

( -2d, d = (i?f - Rl)/(6Rl), R 2 <r< R 3 . 


That solution is displayed in Figure 19. 

Corresponding magnetic flux u(r ) has the following analytical formula 


u(r) 


B 2 — d?i 

2 


d ) r 2 , 


0 < r < Ri, 


— d) r 2 — — — Ri<r<R 2 , 
{ d{R\-r 2 ), R2<r<Rz, 

which can be found as a solution to the following ID problem 


(49) 


_i-i^ = ft „(0)=«(J? 3 ) = 0 (50) 

or r or 

with a piece-wise uniform current density p = 1 for 7?i < r < R 2 and p = 0 otherwise. The coil radii 
are R\ = 0.4, ft 2 = 0.6, and domain size i? 3 = 1.0 for the field shown in Figure 19. 

When the 9FD4U method is used, the computations give zero error (in the sense of accuracy of the 
iterative method) both for the uniform and nonuniform grids. It w r as expected, that the error should 
be of order 0{h 2 ) because of the discontinuity of the right-hand side and nonuniform grid. In fact, 
the error is even better than 0(h 4 ). 

The application of the 5-point scheme 5FD2N (6) for the boundary value problem (50) with the 
same ID solution, but for different f? 3 , provides the following results presented in Figure 20. Here the 
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Figure 19: One-dimensional magnetic field for the infinite solenoid with internal radius 7?i = 0.5 and 
external radius /? 2 = 1.0. Outside conducting wall is assumed to be placed at R 3 — 2.5. 



50 100 200 400 


N 

Figure 20: Relative error of the ID calculation using 5-point scheme ( 6 ) on nonuniform mesh. R 1 — 
0.49, ff 2 - 0.51, = 1, 3, 10. For each N r and i? 3 . the nonuniform mesh was built to minimize the 

relative error. 






relative error is defined as 


b 


max {|«(n) - Mfc(r t -)|} 
0<r,<fi 3 


max 

0<r,<R 3 




( 51 ) 


The plots show that the error is going down as fast as N 2 , where A> is the number of grid points. 
Also, error is less for smaller domain size R 3 due to smaller mesh size. 

Let’s describe the algorithm of generating an adaptive nonuniform mesh for the following ID 
domain in radial direction. The computational domain ft = [0, /? 3 ] was divided into five subdomains: 

ft a = [0, /?i/2], ft 2 = [7?i/2,/?i],ft 3 = [Ru(Ri + fl 2 )/2],ft 4 = [Oh + /? 2 )/2, R&U& = i^Rs}- The 

grid was created to be fine next to singular points: r = 0,r = R\ and r = R 2 . In fti,ft 3 and fts the 
grid step size satisfies the following relation: /*,+ 1 = /ij( 1 + Sh), in ft 2 and ft4: h,+ 1 = h«/( 1 + bh), 
where bh > 0 is a relative increment of the grid-step. This means that in each subdomain ft*, the 
length L(9-k) can be expressed through a number of mesh steps A r £ as following: 


I(ftfc) = K 


,(1 + Sh) + A miw (l + 


t (l + bh) 1 


{(1 + Sh) 1 * - i) /bh. 


The ID adaptive mesh generation has the following three steps. 

T) For given N r number of points and for given relative mesh increment bh, create nonuniform 
mesh in the domain ft. This requires us to find 6 unknowns: the number of grid steps in subdomains 
NINONS, Nliq and the minimal mesh size h m i n . The unknowns must satisfy the following 6 
equations: 

/w((l + Sh) 1 * - 1 ysh = L(ftfc ), k = 1,-5;E N t = ( 52 ) 

k=l 

The following iteration process works efficiently to resolve the above system: 

a) Start with initial guess for A r £ ; 

b) Calculate /i mtn = L(Qs)Sh/((l + bh) N s - 1); 

c) Calculate NL = ln(l + L(Clk)Sh/h m i n )/ ln(l + bh) for k — 1, ...4; 

d) Check: if Y.l=i < N r, then increase A r | by one and return to b). If E/Ui > iV r- then reduce 
A 7 ! by one and return to b) until convergence. 

2) Calculate error b (51) of the numerical solution for the generated mesh, using analytical solution 
«(r) (49). 

3) Optimize mesh by minimizing the error b with respect to the relative mesh increment bh. It 
was found that the function 6(Sh) is a convex smooth function with one minimum (see Figure 21), 
which can be found using the following iteration process: 

a) Start with minimal bh 1 , maximal bh 4 and intermediate bh 3 : bh 1 < bh 3 < bh 4 and calculate 
corresponding errors bi,bz and £4. using st ep 2); 

b) Calculate another intermediate bh 2 = Vb h l bh 3 with corresponding error b 2 \ 

c) If b 2 > S 3 , then assign bh 1 = bh 4 ,bh 4 = bh 2 and return to step b), else assign bh 4 = bh l ,bh 4 = 
bh 3 , bh 3 = bh 2 and return to step b). 

The described nonuniform mesh in radial direction can be efficiently used for solving 2D problems, 
because the 2D solution at the vertical symmetry plane z = 0 behaves similarly to the ID solution 

( 48 ). 

Now let us describe the method of generating an adaptive nonuniform mesh in the axial domain. 
To do this we use a similar approach as used for the radial grid. Consider behavior of the magnetic 
field of a single coil along the middle radius line r = (Ri + R 2 )/ 2. For values of z inside the coil: 
« < . radial component B r goes up from zero to its maximal value at the coil edge. Outside the 

coil, B r goes down, approaching zero at large r ■ It is hard to find a ID magnetostatic problem in 
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Figure 21: Minimization of the numerical error with respect to the relative mesh step increment. 


Cartesian coordinates with a solution having that behavior, but a solution of the ID problem for an 
infinite cylindrical wire with uniform current behaves in a very similar way: 


B(z) 



0 < z < Z 2 , 
Z 2 < z < Z 3 , 


(53) 


That solution is displayed in Figure 22. 

Corresponding magnetic potential .4( r) has the following analytical formula 


A(x) 


Z 2 - 


+ 


7 2 

Z2.lv 


Z J., n (h 



, 0 < z < Z 2 , 
Z 2 <z< Z3, 


(54) 


which can be found as a solution to the following ID Poisson problem 

(S5) 

with a piece-wise uniform charge density function p = 1 for 0 < 2 < Z 2 and p — 0 otherwise. 

In formulas (53, 54, 55) we used variable z for the radial coordinate because the formulas are used 
to generate an adaptive ~-mesh as described below. 

The computational domain 0 = [0,Z3] was divided into two subdomains: Di = [0. Z 2 ],fl 2 = 
\Z 2 . Z3]. The grid was created to be fine next to the only singular point: 2 = Z2. I11 flj the grid step 
size satisfies to the following relation: /?,-+ 1 = /?;/(! + M), in Q 2 \ /) t + 1 = /#»( 1 + fth), where 6 h > 0 is 
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Figure 22: One-dimensional magnetic field for the infinite wire with radius Z 2 = 0.5. 

a relative increment of the grid-step. It means that in each subdomain 0^, the length L(Q k ) can be 
expressed through a number of mesh steps JV£ as following: 

L(Ofc) = h mm + h min ( 1 + 6h) + It mini 1 + 6h) 2 + ...h min { 1 + bhft- 1 = h mi n ((1 + 6h) N l - l) /6h. 

The ID adaptive mesh generation has the following three steps. 

1) For given N r number of points and for given relative mesh increment bh, create a nonuniform 
mesh in the domain fl. This requires us to find 3 unknowns: the number of grid steps in subdomains 
A T j\ AT and the minimal mesh size h min . The unknowns must satisfy the following 3 equations: 

h min (( i + shfk - i; ysh = L(n k ), k = 1,2; at + at = AV; ( 56) 

The following iteration process works efficiently to resolve the above system: 

a) Start with initial guess for AT; 

b) Calculate h m i n — Z(Q 2 )^///((l + 6h) N * - 1); 

c) Calculate A T f = ln(l + L(Qi)6h/h min )/ ln(l + Sh)\ 

d) Check: if V[ + AT < Ay, then increase AT by one and return to b). If jVf + AT > A n then reduce 
AT by one and return to b) until convergence. 

2) Calculate error b (51) of the numerical solution for the generated mesh, using the analytical 
solution u(r) (5-1). 

3) Optimize the mesh by minimizing the error b with respect to the relative mesh increment bh . 
The minimization is done using the same method as for generating the adaptive r-rnesh, described 
above. 
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Appendix 5: Calculation of the magnetic field on the symmetry line 


Accurate calculation of the magnetic field B z on the symmetry line requires use of both t he fine grid in 
the radial direction and the high-order finite difference approximation. Taking into consideration the 
zero conditions for magnetic flux: u(0,~) = 0 and the magnetic potential: .4(0,;) = 0, the following 
analytical formula can be used for computing the axial magnetic field: 




The formula can be approximated be the following second-order finite difference expression: 


(*;: 


or the fourth-order finite difference approximation: 


B h z { 0,~i) 


+ h 


h\ + hi 


Ml ,j 


(_ji_ 

\(*i + i< r 2 : 


“2 ,j 


+ 0(h 4 ). 


(57) 


(58) 


These approximations are derived from the Taylor expansion of u(r, z) near the axis under the condition 
u.\ j — 0, as well as the smoothness and symmetry property of the solution. Let us note that for a test 
function u — r 2 + or 4 , the exact value is B* — 2 and its approximations B % of the second and fourth 
orders equal 2(1 + ah 2 ) and 2 respectively. 


Appendix 6: Convergence in the boundary condition calculation 


The numerical error of the magnetic field depends on several factors: the distance H ~ min{Z 3 — 
R 3 —R 2 } of the domain boundary from the coil, the number m of quadrature points in the numerical 
integration (43), the number n of the roundoff in the series (42) and the characteristic mesh-step 
(R 3 “ # 2 )/™- We consider the influence of these parameters on the example of the rectangular coil 
cl. 

We use three variants of the domain boundary at a different distance from the coil. These bound- 
aries are defined by the following computational domains (Q = {0 < r < R 3 , 0 < z < Z 3 }): 


fti : Rs = ^3 = 2, : R 3 = Z 3 = 4, ft 3 : R 3 = ^3 - 8. 


Tfere we take into account that, due to the symmetry property of the field for the single coil, the 
Neumann boundary condition du/8z = 0 is used at z - 0 and the computational domain includes the 
half of the coil for z > 0. We need to distinguish t hese domains because the value of .s 2 decreases when 
the distance H increases. So in the series (42) we can take a smaller number of terms n and decrease 
the computational costs. Moreover, the underintegral function in (43) is smoother in this case. So it 
is possible to use a smaller number m in the quadrature formulas. 

In Table 3 we present the relative errors 


6 = max 
3 



u a (lh,Zj) - zj) 

U a (R 3,Zj) 


of the boundary values for different integers m, distances i? 3 , and sufficiently big n = n. Errors of 
numerical integration are denoted by 6 r and for the rectangular quadrature formula (43) and the 
Simpson’s quadrature formula (41) respectively. Actually, we define 
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whore m is the minimal value of m which satisfies the condition 


max 

i 




< C = 10 


-13 


It means that, in some sense, is an ""almost exact numerical solution at the boundary. 

Here m = 200, 77 = 32 were taken. 


m 

1 

2 

4 

8 

16 

ST i ; 

S s 

1.M0 -4 

8.2-10 -9 

5 . 7 - 1 0 — 7 

CO 

o 

1 

oc 

2.3-10^ 

: 

S y 

0.045 

0.011 

0.0028 

0.00069 

0.00017 

ST 2 • 

Ss 

2 . 7 - 1 0 — 5 

1.7-i0“ 6 

i . l - io — 7 

6.7-10 -9 

4.2-10 -10 

9. 2 : 

Sr 

0.038 

0.0095 

0.0024 

0.00060 

0.00015 

9 3 -. 

S s 

5.6-10 -6 

3.510 -7 

2.2- 10 -8 

1.4-10 -9 

8.6-10 -11 

9s: 

Sr 

0.036 

0.0091 

0.0023 

0.00057 

0.00014 


Table 3. Errors of numerical integration for coil “cl” 


The table shows that the error is of order 0(m 2 ) for the central rectangle and 0{m 4 ) for the 
Simpson's quadratures and the last ones give an essentially smaller error. For a big distance II the 
quadratures are more precise. Let us note that the calculations using Gauss quadratures provide 
considerably higher efficiency. Its error for m — 7 satisfies the inequality S g < 10 14 for all three types 
of the computational domain. 

Table 4 includes the errors 

(n 3 . =J )-u^(P 3 .z 3 ) 1 

J 

of the boundary values for a different number of series terms in (42). The quadrature used is Simpson’s 
one. Here u a {R 3 .Zj ) is the “e-asymptotically” exact numerical solution with m - 200, n - 64. 


n 

2 

4 

8 

16 

32 

^3 

0.44 

0.19 

0.067 

0.22 

0.045 

0.0056 

0.072 
0.0034 
5.8- 10 -5 

0.011 
3.5- 10- 5 
1.2- 10 -8 

5.9-1 O' 4 
8.4-1 0' 9 
2.0-10 -15 


Table 4. Errors of computing the elliptic integrals for coil “cl” 
Note that n = 32 is not enough for if high accuracy is required. 
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